#include <iostream>
#include <iomanip>
#include <string>
#include <numeric>
#include <map>
#include <ginac/ginac.h>
#include <boost/lexical_cast.hpp>
#include <boost/utility.hpp>
#include <boost/foreach.hpp>
#include <boost/tuple/tuple.hpp>

#include <boost/fusion/container/map.hpp>
#include <boost/fusion/include/map.hpp>
#include <boost/fusion/container/map/map_fwd.hpp>
#include <boost/fusion/include/map_fwd.hpp>
#include <boost/fusion/sequence/intrinsic/at_key.hpp>
#include <boost/fusion/include/at_key.hpp>
#include <boost/fusion/container/generation/make_map.hpp>
#include <boost/fusion/include/make_map.hpp>



#include <lemon/lp.h> // linear programming solver
#include <boost/random.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/math/special_functions/modf.hpp>



#include <boost/icl/interval.hpp>
#include <boost/icl/discrete_interval.hpp>
#include <boost/icl/continuous_interval.hpp>
#include <boost/icl/right_open_interval.hpp>
#include <boost/icl/left_open_interval.hpp>
#include <boost/icl/closed_interval.hpp>
#include <boost/icl/open_interval.hpp>
#include <boost/icl/functors.hpp>


using namespace GiNaC;
using namespace boost;
using std::cout;
using std::setw;
using std::left;
using std::endl;
using std::string;
using std::distance;
using std::make_pair;
using std::pair;
using std::accumulate;

/// Factory for symbol creation
const symbol & get_symbol(const string & s)
{
  static std::map<string, symbol> directory;
  std::map<string, symbol>::iterator i = directory.find(s);
  if (i != directory.end())
    return i->second;
  else
    return directory.insert(make_pair(s, symbol(s))).first->second;
}
/// Enum for easy numeration of UFx_lst_tuple entries
enum {Uarg,Farg,x_lstarg};
/// Tuple for storing U and F polynoms and list of Feynman parameters x_lst
typedef tuple<ex,ex,lst>  UFx_lst_tuple;
namespace UFX
{
  struct U;
  struct F;
  struct xlst;
}

typedef fusion::map<fusion::pair< UFX::U, ex>, fusion::pair< UFX::F, ex>, fusion::pair< UFX::xlst,lst> > UFXmap;








std::pair<ex,double>  hyper_cube(lst pole_list,lst w_list)
{
  try{
    using namespace lemon;
    Lp lp;
    exhashmap<LpBase::Col> col_map;
    for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      {
        col_map[*wi] = lp.addCol();
      }
    for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
      {
        ex tmp_expr = *pit;
        // cout<<"Expr: "<<*pit<<" subexprs:"<<endl;
        Lp::Expr constr_expr;
        for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
          {
            // cout<<*wi<<"  coeff "<<tmp_expr.coeff(*wi)<<endl;
            ex wi_coeff = tmp_expr.coeff(*wi);
            tmp_expr-=(*wi)*wi_coeff;
            if(is_a<numeric>(wi_coeff))
              {
                constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
                //                cout<<ex_to<numeric>(wi_coeff).to_double()<<endl; 
              }
            else throw "Non numeric coefficient in pole term. ";
          }
        //constr_expr+=
        //cout<<"Ostatok "<<tmp_expr<<endl;
        if(is_a<numeric>(tmp_expr))
          lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
        else throw string("Lower bound is not a numeric");
      }

    double l_bound,r_bound;
    cout<<"Hyper cube"<<endl;
    //   for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      
        lp.min();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<*w_list.begin()<<" = ("<<lp.primal()<<",";
            l_bound = lp.primal();
          }
        else throw  string("Optimal solution not found.");
        lp.max();
        lp.obj(col_map[*w_list.begin()]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          {
            cout<<lp.primal()<<")"<<endl;
            r_bound = lp.primal();
          }
        else throw  string("Optimal solution not found.");

        boost::mt19937 rng;       
        boost::uniform_real<> six(l_bound,r_bound);      // distribution that maps to 1..6
        // see random number distributions
        boost::variate_generator<boost::mt19937&, boost::uniform_real<> >
          die(rng, six);             // glues randomness with mapping
        rng.seed(static_cast<unsigned char>(std::time(0)));

        
        double x   = die();                      // simulate rolling a die

        cout<<"Point "<<die()<<"   "<<x<<endl;
        return std::make_pair(*w_list.begin(),x);
        
  }
  catch(string& what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
}


exmap start_point(lst pole_list,lst w_list)
{
  exmap subs_map;
  for(lst::const_iterator wi = w_list.begin();wi != w_list.end();++wi) 
    {
      lst tmp_pole_list;
      lst tmp_w_list;
      for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
        if(!((*pit).subs(subs_map)>0))
          tmp_pole_list.append((*pit).subs(subs_map));
      cout<<"tmp_pole_list"<<tmp_pole_list<<endl;
      for(lst::const_iterator wi2 = wi;wi2 != w_list.end();++wi2) 
        tmp_w_list.append(*wi2);
      std::pair<ex,double> ret_pair = hyper_cube(tmp_pole_list,tmp_w_list);
      subs_map[ret_pair.first] = ret_pair.second;
      cout<<ret_pair.first<<" "<<ret_pair.second<<endl;
    }
  cout<<subs_map<<endl;
  return subs_map;
}


class MBintegral
{
  ex no_gamma_part;
  lst gamma_poles;
  lst w_lst;
  exmap eps_w_current;
  exmap w_current;
  relational eps_current;
public:
  typedef std::list<double>::iterator gamma_iterator;
  typedef lst::const_iterator pole_iterator;
  std::pair<pole_iterator,pole_iterator> gamma_args()
  {
    return std::make_pair(gamma_poles.begin(),gamma_poles.end());
  }
  lst get_pole_lst()
  {
    return gamma_poles;
  }
  lst get_w_lst()
  {
    return w_lst;
  }
  exmap new_point()
  {
    lst var_list(w_lst);
    var_list.append(get_symbol("eps"));
    eps_w_current=start_point(get_pole_lst(),var_list);
    for(lst::const_iterator it = w_lst.begin();it!=w_lst.end();++it)
      w_current[*it] = it->subs(eps_w_current);
    eps_current = (get_symbol("eps")==eps_w_current[get_symbol("eps")]);
    return eps_w_current;
  }

  exmap get_point()
  {
    return eps_w_current;
  }
  exmap get_w()
  {
    return w_current;
  }
  relational get_eps()
  {
    return eps_current;
  }
  MBintegral res(relational w_relation,ex pole,relational new_eps)
  {
    lst new_w_lst;
    for(lst::const_iterator it = w_lst.begin();it!=w_lst.end();++it)
      if(*it != w_relation.lhs()) new_w_lst.append(*it);
    //    std::remove_copy(w_lst.begin(),w_lst.end(),new_w_lst.begin(),w_relation.lhs());
    cout<<"removed "<<w_relation.lhs()<<"  in list "<<new_w_lst<<endl;
    lst new_gamma_pole_list;
    for(lst::const_iterator it = gamma_poles.begin();it!=gamma_poles.end();++it)
      if(*it != pole) new_gamma_pole_list.append(it->subs(w_relation));
    cout<<"removed pole list  "<<pole<<"  in list "<<new_gamma_pole_list<<endl;
    // !!! IMPORTANT!!! no 2*pi*i multiplication and no sign multiplication 
    ex new_no_gamma_part = (no_gamma_part.subs(tgamma(pole)==pow(-1,pole.subs(w_relation))/factorial(pole.subs(w_relation)))).subs(w_relation);
    // new_no_gamma_part  = pow(-1,pole.subs(w_relation))/factorial(pole.subs(w_relation))*no_gamma_part.subs(w_relation);
    cout<< new_no_gamma_part<<endl;
    exmap new_w_current(w_current);
    cout<<" Not modif:  "<<new_w_current<<endl;
    new_w_current.erase(w_relation.lhs());
    cout<<" modif:  "<<new_w_current<<endl;
    MBintegral resINT(new_w_lst,new_gamma_pole_list,new_no_gamma_part,new_w_current,new_eps);
    return resINT;
  }
   ex get_expr()
  {
    return no_gamma_part; 
  }

  ex get_gamma_expr()
  {
    ex tmpex = no_gamma_part;
    for(pole_iterator it = gamma_poles.begin();it!=gamma_poles.end();++it)
      {
        tmpex*=tgamma(*it);
      }
    return tmpex;
  }
 ex operator*=(ex to_mul)
  {
    no_gamma_part*=to_mul;
    return no_gamma_part;
  }
  // Constructor for Residue
  MBintegral(lst w_lst_in,lst pole_lst_in,ex no_gamma_part_in,exmap w_current_in,relational eps_in):w_lst(w_lst_in),gamma_poles(pole_lst_in),no_gamma_part(no_gamma_part_in),w_current(w_current_in),eps_current(eps_in)
  {
 
  }

  MBintegral(UFXmap ufx_in,lst nu,numeric l) // lst nu is a list of powers of propagators and l is a number of loops
  {
    ex N = accumulate(nu.begin(),nu.end(),ex(0));
    cout<<"N(nu)= "<<N<<endl;
    lst  x_lst(fusion::at_key<UFX::xlst>(ufx_in));
    ex U = fusion::at_key<UFX::U>(ufx_in).collect(x_lst,true);
    ex F = fusion::at_key<UFX::F>(ufx_in).collect(x_lst,true);///< distributed polynom factorizing X(i)X(j)*(...)
    
    cout<<setw(30)<<std::internal<<"+++INTEGRAL PARAMETERS+++"<<endl;
    cout<<setw(30)<<std::left<<"** Number of loops L=   "<<std::internal<<l<<endl;
    cout<<setw(30)<<std::left<<"** Summ of powers N=   "<<N<<endl;
    //cout<<setw(35)<<std::left<<"** Dimension D=   "<<D.subs(D_subs)<<endl;
    cout<<setw(35)<<std::left<<"** U =   "<<U<<endl;
    cout<<setw(35)<<std::left<<"** F =   "<<F<<endl;
    
    // assuming 1/(U^a*F^b)
    ex U_pow = (-N+(l+1)*(2-get_symbol("eps")));
    cout<<"U_pow"<<U_pow<<endl;
    ex F_pow = (N-l*(2-get_symbol("eps")));
    //lst gamma_lst; //Gamma with poles left of contour
    //    lst w_lst;
    //  ex gamma_right; // and right of contour
    exmap x_power_map; // map of X(j) powers
    //    ex coeff = 1;                             // numerical coeeficient independent of X(j)
    ex coeff = pow(exp(get_symbol("eps")*Euler)*pow(get_symbol("mu"),2*get_symbol("eps")),l)*tgamma(F_pow);
    // important if power = 0????
    for(lst::const_iterator nui = nu.begin();nui!=nu.end();++nui)
      coeff/=tgamma(*nui);
    coeff/=tgamma(F_pow);
    //    ex out_ex = 1;///pow(2*Pi*I,U.nops()-1)/tgamma(al_pow)/pow(2*Pi*I,F.nops()-1)/tgamma(al_pow); //need review
    //working with F-term \Gamma(\nu-L*D/2) contractedx
    ex w_sum = 0;  //F-term generates only integrations in W
    for(const_iterator it = F.begin();it!=F.end();++it)
      {
        int w_index = distance(F.begin(),it);
        ex x_power;
        if(F.end()==boost::next(it)) 
          {
            coeff*=tgamma(F_pow+w_sum);
            x_power = -F_pow-w_sum;
            gamma_poles.append(F_pow+w_sum);
            cout<<"end achieved"<<endl;
          }
        else
          {
            string str = "w_"+boost::lexical_cast<string>(w_index);
            //symbol w(str);
            w_lst.append(get_symbol(str));
            coeff*=tgamma(-get_symbol(str));
            x_power = get_symbol(str);
            gamma_poles.append(-get_symbol(str));
            w_sum+=get_symbol(str);
            cout<<"ok run"<<endl;
          }
        // filling map of X(j) powers
        ex tmp_expr = (*it);
        for(lst::const_iterator it1=x_lst.begin();it1!=x_lst.end();++it1)
          {
            if(tmp_expr.has(*it1))
              {
                x_power_map[*it1]+=(tmp_expr.degree(*it1)*x_power);
                tmp_expr = tmp_expr.subs((*it1)==1);
              }
          }
        // coeff *=pow(tmp_expr,x_power);
      }
    
    //working with U-term, only if (U_pow>0)
    
    
    if(U_pow.subs(get_symbol("eps")==0) > 0)
      {
        cout<<">>> Working with MB for U term, U_power = "<< U_pow<<endl;  
        ex v_sum = 0;  //F-term generates only integrations in V
        for(const_iterator it = U.begin();it!=U.end();++it)
          {
            int v_index = distance(U.begin(),it);
            ex x_power;
            if(U.end()==boost::next(it)) 
              {
                //          out_ex*=pow(*it,-al_pow-w_sum)*tgamma(al_pow+w_sum);
                x_power = -U_pow - v_sum;
                gamma_poles.append(U_pow+v_sum);
                cout<<"end achieved"<<endl;
              }
            else
              {
                string str = "v_"+boost::lexical_cast<string>(v_index);
                //symbol v(str);
                w_lst.append(get_symbol(str));
                //  out_ex*=pow(*it,w)*tgamma(-w);
                x_power = get_symbol(str);
                gamma_poles.append(-get_symbol(str));
                w_sum+=get_symbol(str);
                cout<<"ok run"<<endl;
              }
            // filling map of X(j) powers
            ex tmp_expr = (*it);
            for(lst::const_iterator it1=x_lst.begin();it1!=x_lst.end();++it1)
              {
                if(tmp_expr.has(*it1))
                  {
                  x_power_map[*it1]+=(tmp_expr.degree(*it1)*x_power);
                  tmp_expr = tmp_expr.subs((*it1)==1);
                }
              }
            coeff *=pow(tmp_expr,x_power);
          }
      }
    else
      {
        ex U_term = pow(U,U_pow.subs(get_symbol("eps")==0));
        cout<< "U_TERM  "<<U_term<<endl;
      }
    
    //cout<<x_power_map<<endl;
    cout<<"Gammas after MB: "<<endl<<gamma_poles<<endl;
    cout<<"X powers list:"<<endl<<x_power_map<<endl;
    // applying X integration
    ex gamma_den = 0; // gamma in denominator
    for(exmap::const_iterator mi = x_power_map.begin();mi!=x_power_map.end();++mi)
      {
        cout<<(*mi).first<<" "<<(*mi).second<<endl;
        gamma_poles.append((*mi).second+1);
        coeff*=tgamma((*mi).second+1);
        gamma_den+=(*mi).second+1;
      }
    coeff/=tgamma(gamma_den);
    cout<<"New gamma list:"<<endl<<gamma_poles<<endl;
    no_gamma_part = coeff;
  }
};


lst has_w(ex gamma_arg,lst w_list)
{
  lst out_list;
  for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
    if(gamma_arg.has(*wi))out_list.append(*wi);
  return out_list;
}


void hit_pole(lst pole_list,lst w_list,exmap subs_map)
{
  lst C_list;
  exmap C_subs_map;
  for(lst::const_iterator cit = w_list.begin();cit!=boost::prior(w_list.end()); ++cit)
    {
      C_list.append(*cit);
      C_subs_map[*cit] = subs_map[*cit];
    }

  ex eps0 = *w_list.rbegin();


  cout<<eps0.subs(subs_map)<<"  "<<boost::math::iround(ex_to<numeric>(eps0.subs(subs_map)).to_double())<<endl;
  cout<<boost::math::iround(3.7)<<" 3.7 "<<boost::math::iround(-3.7)<< " -3.7"<<endl;
  cout<<"poles closest to C "<<endl;
  for(lst::const_iterator pit = C_list.begin();pit!=C_list.end();++pit)
    {
      int ipart = boost::math::iround(ex_to<numeric>((*pit).subs(subs_map)).to_double());
      double fpart = ex_to<numeric>((*pit).subs(subs_map)).to_double();
      double delta = std::abs(fpart - ipart);

      cout<< "("<<*pit<<","<<ipart<<")   "<<" float: "<<fpart<<" ipart: "<<ipart<<" delta: "<<delta<<endl;
      

    }
  cout<<endl;
  ex eps_max = 0;

  for(lst::const_iterator ui = pole_list.begin();ui!=pole_list.end();++ui)
    {
      lst U_w_list(has_w(*ui,C_list));
      if((U_w_list.nops()>0)&&(ui->has(eps0)))
        {
          cout<<U_w_list<<endl;
          std::vector<double> delta_vector;
          double min_delta = 1000;
          int U_min;
          for(lst::const_iterator wi=U_w_list.begin();wi!=U_w_list.end();++wi)
            if(std::abs(ex_to<numeric>((*wi).subs(subs_map)).to_double() - boost::math::iround(ex_to<numeric>((*wi).subs(subs_map)).to_double()))<min_delta) 
              {
                min_delta = std::abs(ex_to<numeric>((*wi).subs(subs_map)).to_double() - boost::math::iround(ex_to<numeric>((*wi).subs(subs_map)).to_double()));
                U_min = boost::math::iround(ex_to<numeric>((*wi).subs(subs_map)).to_double());
              }
          cout<<"n(i)="<<U_min<<endl;
          ex tmp_eps = lsolve((*ui).subs(C_subs_map) == U_min,eps0 );
          if(abs(tmp_eps) > abs(eps_max))eps_max = tmp_eps;
          cout<<"EPS1 = "<< tmp_eps<<endl;
          
        }
    }

  cout<<"EPS(1) = "<<eps_max<<endl;
  lst u_n_eqs;
  for(lst::const_iterator ui = pole_list.begin();ui!=pole_list.end();++ui)
    {
      // pole value in new EPS point:
      ex U_pole = (ui->subs(C_subs_map)).subs(eps0 == eps_max);
      cout<<*ui<<"         "<<U_pole<<endl;
      if(U_pole==0) u_n_eqs.append(*ui ==0);
      else if(U_pole==-1) u_n_eqs.append(*ui ==-1);
      else if(U_pole==-2) u_n_eqs.append(*ui ==-2);
      else if(U_pole==-3) u_n_eqs.append(*ui ==-3);
      else if(U_pole==-4) u_n_eqs.append(*ui ==-4);
      else if(U_pole==-5) u_n_eqs.append(*ui ==-5);
    }
  // IF P{...}!=0
  if(u_n_eqs.nops()>0)
    {
      lst w_in_eq;
      for(lst::const_iterator eqi = u_n_eqs.begin();eqi!=u_n_eqs.end();++eqi)
        {
          lst wlst = has_w(eqi->lhs(),C_list);
          cout<<wlst<<endl;
          
          for(lst::const_iterator wi = wlst.begin();wi!=wlst.end();++wi)
          w_in_eq.append(*wi);
        }
      w_in_eq.unique();
      ex w_sub_solve = lsolve(u_n_eqs,w_in_eq);
for(lst::const_iterator ui = pole_list.begin();ui!=pole_list.end();++ui)
      cout<<ui->subs(w_sub_solve).subs(C_subs_map).subs(eps0 == eps_max)<<endl;
    }
}

void gamma_lp(lst pole_list,lst w_list)
{
  try{
    using namespace lemon;
    Lp lp;
    exhashmap<LpBase::Col> col_map;
    for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      {
        col_map[*wi] = lp.addCol();
      }
    for(lst::const_iterator pit = pole_list.begin();pit!= pole_list.end();++pit)
      {
        ex tmp_expr = *pit;
        // cout<<"Expr: "<<*pit<<" subexprs:"<<endl;
        Lp::Expr constr_expr;
        for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
          {
            // cout<<*wi<<"  coeff "<<tmp_expr.coeff(*wi)<<endl;
            ex wi_coeff = tmp_expr.coeff(*wi);
            tmp_expr-=(*wi)*wi_coeff;
            if(is_a<numeric>(wi_coeff))
              {
                constr_expr+=ex_to<numeric>(wi_coeff).to_double()*col_map[*wi];
                //                cout<<ex_to<numeric>(wi_coeff).to_double()<<endl; 
              }
            else throw "Non numeric coefficient in pole term. ";
          }
        //constr_expr+=
        //cout<<"Ostatok "<<tmp_expr<<endl;
        if(is_a<numeric>(tmp_expr))
          lp.addRow(-ex_to<numeric>(tmp_expr).to_double(),constr_expr,Lp::INF);
        else throw "Lower bound is not a numeric";
      }

    cout<<"Hyper cube"<<endl;
    for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
      {
        lp.max();
        lp.obj(col_map[*wi]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          cout<<*wi<<" = ("<<lp.primal()<<",";
        lp.min();
        lp.obj(col_map[*wi]);
        // Solve the problem using the underlying LP solver
        lp.solve();
        if (lp.primalType() == Lp::OPTIMAL) 
          cout<<lp.primal()<<")"<<endl;
      }

    cout<<"Maximal eps value"<<endl;
  // Specify the objective function
  lp.max();
  lp.obj(col_map[*w_list.rbegin()]);
  
  // Solve the problem using the underlying LP solver
  lp.solve();

  // Print the results
  if (lp.primalType() == Lp::OPTIMAL) 
    {
std::cout << "Objective function value: " << lp.primal() << std::endl;
      for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
        {
            std::cout <<*wi<< " = " << lp.primal(col_map[*wi]) << std::endl;
        }
    } 
  else throw  "Optimal solution not found.";

    cout<<"Minimal eps value"<<endl;
  // Specify the objective function
  lp.min();
  lp.obj(col_map[*w_list.rbegin()]);
  
  // Solve the problem using the underlying LP solver
  lp.solve();

  // Print the results
  if (lp.primalType() == Lp::OPTIMAL) 
    {
std::cout << "Objective function value: " << lp.primal() << std::endl;
      for(lst::const_iterator wi = w_list.begin();wi!=w_list.end();++wi)
        {
            std::cout <<*wi<< " = " << lp.primal(col_map[*wi]) << std::endl;
        }
    } 
  else throw  "Optimal solution not found.";
  

  }
  catch(string what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
}


UFXmap UF(lst k_lst,lst p_lst,lst subs_lst)
{
  /// constructing expr Xi*Pi
  ex fprop = 0;
  ex sDtmp;
  string str;
  lst xp;
  for(lst::const_iterator it = p_lst.begin();it!=p_lst.end();++it)
    {
      str = "x_" + boost::lexical_cast<string>(std::distance(p_lst.begin(),it));
      xp.append(get_symbol(str)); // storing list of X identities
      fprop += get_symbol(str) * (*it); 
    }
  fprop = fprop.expand();
  sDtmp = fprop;
  cout<<fprop<<endl;
  matrix M(k_lst.nops(),k_lst.nops());
  matrix Q(k_lst.nops(),1);//column
  GiNaC::numeric half(1,2);
  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    for(lst::const_iterator itc = itr;itc!=k_lst.end();++itc)
      if(itr == itc)
        {
          M(distance(k_lst.begin(),itr),distance(k_lst.begin(),itc)) = -1*fprop.coeff((*itr)*(*itc),1);
          cout<<(*itr)*(*itc)<<"M("<<distance(k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr)*(*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr)*(*itc),1);
        }
      else
        {
          M(distance(k_lst.begin(),itr),distance( k_lst.begin(),itc)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          M(distance(k_lst.begin(),itc),distance(k_lst.begin(),itr)) = -1*half*fprop.coeff((*itr),1).coeff((*itc),1);
          cout<<(*itr)*(*itc)<<"M("<<distance( k_lst.begin(),itr)<<","<<distance( k_lst.begin(),itc)<<") coeff "<<fprop.coeff((*itr),1).coeff((*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*fprop.coeff((*itr),1).coeff((*itc),1);
        }
  cout<<"M: "<<M<<endl;
  sDtmp = sDtmp.expand();
  cout<<"Expr linear on external momentum: "<<sDtmp.expand()<<endl;

  for(lst::const_iterator itr = k_lst.begin();itr!=k_lst.end();++itr)
    {
      Q(distance( k_lst.begin(),itr),0) = half*sDtmp.coeff((*itr),1);
      sDtmp -= (*itr)*sDtmp.coeff((*itr),1);
    }
  cout<<"Q: "<<Q<<endl;
  sDtmp = sDtmp.expand();
  ex minusJ = sDtmp;
  cout<<"-J: "<<minusJ<<endl;
  ex U = M.determinant();
  ex F = expand(M.determinant()*(minusJ+Q.transpose().mul(M.inverse()).mul(Q)(0,0)));
  lst lp;
  F=F.normal();
  cout<<"U= "<<U<<endl<<"F= "<<F<<endl;
  //  cout<<"pol_list "<<lp<<endl;
 U=U.subs(subs_lst);
  F=F.subs(subs_lst);
  return fusion::make_map<UFX::U,UFX::F,UFX::xlst>(U,F,xp);
}


ex D;
//symbol eps("eps");
/// list variant
ex MB_lst(pair<lst,lst> UF_x_lst,lst nu,numeric l,relational D_subs=D==4-2*get_symbol("eps"))
{
  D = D.subs(D_subs);
  cout<<D<<endl;
  ex N = accumulate(nu.begin(),nu.end(),ex(0));
  cout<<"N(nu)= "<<N<<endl;
  lst  x_lst(UF_x_lst.second);
  ex U = UF_x_lst.first.op(0).subs(D_subs);
  ex F = UF_x_lst.first.op(1).subs(D_subs).collect(x_lst,true);///< distributed polynom factorizing X(i)X(j)*(...)

  cout<<setw(30)<<std::internal<<"+++INTEGRAL PARAMETERS+++"<<endl;
  cout<<setw(30)<<std::left<<"** Number of loops L=   "<<std::internal<<l<<endl;
  cout<<setw(30)<<std::left<<"** Summ of powers N=   "<<N<<endl;
  cout<<setw(35)<<std::left<<"** Dimension D=   "<<D.subs(D_subs)<<endl;
  cout<<setw(35)<<std::left<<"** U =   "<<U<<endl;
  cout<<setw(35)<<std::left<<"** F =   "<<F<<endl;

  // assuming 1/(U^a*F^b)
  ex U_pow = (-N+(l+1)*D/2).subs(D_subs);
  cout<<"U_pow"<<U_pow<<endl;
  ex F_pow = (N-l*D/2).subs(D_subs);
  lst gamma_lst; //Gamma with poles left of contour
  lst w_lst;
  //  ex gamma_right; // and right of contour
  exmap x_power_map; // msp of X(j) powers
  ex coeff = 1;                             // numerical coeeficient independent of X(j)
  ex out_ex = 1;///pow(2*Pi*I,U.nops()-1)/tgamma(al_pow)/pow(2*Pi*I,F.nops()-1)/tgamma(al_pow); //need review
  //working with F-term \Gamma(\nu-L*D/2) contracted
  ex w_sum = 0;  //F-term generates only integrations in W
  for(const_iterator it = F.begin();it!=F.end();++it)
    {
      int w_index = distance(F.begin(),it);
      ex x_power;
      if(F.end()==boost::next(it)) 
        {
          //          out_ex*=pow(*it,-al_pow-w_sum)*tgamma(al_pow+w_sum);
          x_power = -F_pow-w_sum;
          gamma_lst.append(F_pow+w_sum);
          cout<<"end achieved"<<endl;
        }
      else
        {
          string str = "w_"+boost::lexical_cast<string>(w_index);
          //symbol w(str);
          w_lst.append(get_symbol(str));
          //  out_ex*=pow(*it,w)*tgamma(-w);
          x_power = get_symbol(str);
          gamma_lst.append(-get_symbol(str));
          w_sum+=get_symbol(str);
          cout<<"ok run"<<endl;
        }
      // filling map of X(j) powers
      ex tmp_expr = (*it);
      for(lst::const_iterator it1=x_lst.begin();it1!=x_lst.end();++it1)
        {
          if(tmp_expr.has(*it1))
            {
              x_power_map[*it1]+=(tmp_expr.degree(*it1)*x_power);
              tmp_expr = tmp_expr.subs((*it1)==1);
            }
        }
      coeff *=pow(tmp_expr,x_power);
    }

  //working with U-term, only if (U_pow>0)


  if(U_pow.subs(get_symbol("eps")==0) > 0)
    {
      cout<<">>> Working with MB for U term, U_power = "<< U_pow<<endl;  
      ex v_sum = 0;  //F-term generates only integrations in V
      for(const_iterator it = U.begin();it!=U.end();++it)
        {
          int v_index = distance(U.begin(),it);
          ex x_power;
          if(U.end()==boost::next(it)) 
            {
              //          out_ex*=pow(*it,-al_pow-w_sum)*tgamma(al_pow+w_sum);
              x_power = -U_pow - v_sum;
              gamma_lst.append(U_pow+v_sum);
              cout<<"end achieved"<<endl;
            }
          else
            {
              string str = "v_"+boost::lexical_cast<string>(v_index);
              //symbol v(str);
              w_lst.append(get_symbol(str));
              //  out_ex*=pow(*it,w)*tgamma(-w);
              x_power = get_symbol(str);
              gamma_lst.append(-get_symbol(str));
              w_sum+=get_symbol(str);
              cout<<"ok run"<<endl;
            }
          // filling map of X(j) powers
          ex tmp_expr = (*it);
          for(lst::const_iterator it1=x_lst.begin();it1!=x_lst.end();++it1)
            {
              if(tmp_expr.has(*it1))
                {
                  x_power_map[*it1]+=(tmp_expr.degree(*it1)*x_power);
                  tmp_expr = tmp_expr.subs((*it1)==1);
                }
            }
          coeff *=pow(tmp_expr,x_power);
        }
    }


//cout<<x_power_map<<endl;
cout<<"Gammas after MB: "<<endl<<gamma_lst<<endl;
 cout<<"X powers list:"<<endl<<x_power_map<<endl;
// applying X integration
 ex gamma_den = 0; // gamma in denominator
 for(exmap::const_iterator mi = x_power_map.begin();mi!=x_power_map.end();++mi)
   {
     cout<<(*mi).first<<" "<<(*mi).second<<endl;
     gamma_lst.append((*mi).second+1);
     gamma_den+=(*mi).second+1;
   }
 cout<<"New gamma list:"<<endl<<gamma_lst<<endl;
 lst lp_lst(w_lst);
 lp_lst.append(get_symbol("eps"));

 gamma_lp(gamma_lst,lp_lst);
 //cout<<hyper_cube(gamma_lst,lp_lst).rhs();
 exmap w0_map(start_point(gamma_lst,lp_lst));
 cout<<"gamma den:"<<gamma_lst.subs(w0_map)<<endl;
 hit_pole(gamma_lst,lp_lst,w0_map);

return out_ex*coeff;
}

typedef std::list<MBintegral> MBlst;

void MBcontinue(MBintegral rootint,double eps0 = 0)
{
  MBlst O(1,rootint);
  MBlst C;
  while(O.size()>0)
    {
      MBlst R;
      for(MBlst::iterator it = O.begin();it!=O.end();++it)
        {
          C.push_back(*it);//need review, multiple entries C=C U I
          MBintegral::pole_iterator pit,pit_end;
          ex eps_i = get_symbol("eps");
          eps_i = eps_i.subs(it->get_eps());
          cout<<"eps_i = "<<eps_i<<endl;
          for(tie(pit,pit_end) = it->gamma_args();pit!=pit_end;++pit) // loop over gamma arguments
            {
              cout<<"F(eps_i) "<<pit->subs(it->get_w()).subs(it->get_eps())<<"F(eps=0) "<<pit->subs(it->get_w()).subs(get_symbol("eps")==eps0)<<"   min  "<<std::min(pit->subs(it->get_w()).subs(it->get_eps()),pit->subs(it->get_w()).subs(get_symbol("eps")==eps0))<<endl;
              ex F_eps0 = pit->subs(it->get_w()).subs(get_symbol("eps")==eps0);
              ex F_epsi = pit->subs(it->get_w()).subs(it->get_eps());
              for(int n = 0;n>std::min(F_eps0,F_epsi);n--)
                {
                  ex eps_prime = lsolve(pit->subs(it->get_w()) ==n,get_symbol("eps") );
                  cout<<"F= "<<*pit<<endl;
                  cout<<"eps_i: "<<lsolve(pit->subs(it->get_w()) ==n,get_symbol("eps") )<<endl;
                  cout<<" Poles of Gamma on eps_i: "<<it->get_pole_lst().subs(it->get_w()).subs(get_symbol("eps")==eps_prime)<<endl;
                  lst w_in_F  = has_w(*pit,it->get_w_lst());
                  if(w_in_F.nops()>0)
                    {
                      cout<<lsolve(*pit==n,w_in_F.op(0))<<endl;
                      cout<<"sign(z) = "<<csgn(pit->coeff(w_in_F.op(0)))<<"     sign(F_i-F_0) = "<<csgn(F_epsi-F_eps0)<<endl;
                      
                      //MBintegral newi(it->get_w_lst(),it->get_pole_lst(),it->get_expr(),it->get_w(),it->get_eps());
                      //      cout<<"NEW INT: "<< newi.get_expr()<<endl;
                      MBintegral res_int(it->res(w_in_F.op(0)==lsolve(*pit==n,w_in_F.op(0)),*pit,get_symbol("eps")==eps_prime));
                      cout<<"Storing RES_INT wit eps_prime = " <<eps_prime<<"  "<<res_int.get_eps().rhs()<<endl;
                      res_int*=(2*Pi*I*csgn(pit->coeff(w_in_F.op(0)))*csgn(F_epsi-F_eps0));
                      cout<<"RES EXPR:  "<<res_int.get_expr()<<endl;
                      R.push_back(res_int);
                    }
                  else cout<<endl<<endl<<"EEEEEERRRRRRRROOOORR: no W dependence in pole"<<endl<<endl;
                }
            }
        }
      O = R;
    }
  cout<<"Continue get "<<C.size()<<" integrals"<<endl;
  for(MBlst::iterator it = C.begin();it!= C.end();++it)
    {
      cout<<(it->get_pole_lst().subs(it->get_w())).subs(get_symbol("eps")==0)<<endl;
      cout<<endl <<"series:  " << it->get_gamma_expr().series(get_symbol("eps")==0,0)<<endl<<endl;
    }
}

int main()
{
  symbol k("k"),q("q"),p("p"),p1("p1"),p2("p2"),p3("p3"),ms("m^2");
    UFXmap l45 = UF(lst(k),lst(pow(k,2),pow(k+p1,2),pow(k+p1+p2,2),pow(k+p1+p2+p3,2)),lst(pow(p1,2)==0,pow(p2,2)==0));
  MBintegral root_int(l45,lst(1,1,1,1),1);

  //UFXmap l45 = UF(lst(k),lst(ms-pow(k,2),-pow(p1-k,2)),lst(ms==1));
  //MBintegral root_int(l45,lst(1,1),1);

  cout<<" new point "<<endl<<root_int.new_point()<<endl;
  cout<<" saved point "<<endl<<root_int.get_point()<<endl;
  MBcontinue(root_int);
  //cout<<MB_lst(l45,lst(1,1,1,1),1).expand()<<endl;
  return 0;
}





